/*
 * rootGSL.cpp
 *
 *  Created on: 11-09-2012
 *      Author: mith
 */

#include "rootGSL.h"

int RootGSL::Solve(Variables &InitV) {
	const gsl_multiroot_fsolver_type *T;
	gsl_multiroot_fsolver *s;

	unsigned N_variables = InitV.getNtabD();
	int status;
	size_t iter = 0;

	gsl_multiroot_function f = {fun, N_variables, P};

	gsl_vector *x = gsl_vector_alloc (N_variables);
	for(unsigned i=0; i<N_variables; i++) { gsl_vector_set(x, i, InitV[i]); }

	T = gsl_multiroot_fsolver_hybrids;
	s = gsl_multiroot_fsolver_alloc (T, N_variables);
	gsl_multiroot_fsolver_set (s, &f, x);

	//if(DEBUG2) { P->print_label(stdout); InitV.print_label(stdout); print_state (iter, s, N_variables); }

	do {
		iter++;
	    status = gsl_multiroot_fsolver_iterate (s);
		if(DEBUG) { print_state (iter, s, N_variables); }

	    if (status) { break; }  // check if solver is stuck
	    status = gsl_multiroot_test_residual (s->f, PRECISION);
	}
	while (status == GSL_CONTINUE && iter < MAX_N_ITERATIONS);

	if(DEBUG) { printf ("status = %s\n", gsl_strerror (status)); }

	R->set(s);

	gsl_multiroot_fsolver_free (s);
	gsl_vector_free (x);

	return status;
}

void RootGSL::print_state (size_t iter, gsl_multiroot_fsolver * s, unsigned N_variables)
{
	printf("%3u", (unsigned int)iter);
	for(unsigned i=0; i<N_variables; i++) {
		printf(" %15.6e (%3d)", gsl_vector_get (s->x, i), logabs10(gsl_vector_get (s->f, i)));
	}
	printf("\n");
}

void RootGSL::print(FILE *file) const {
	R->print(file);
}

void RootGSL::print_verboten(FILE *file) const {
	R->print_verboten(file);
}
